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We study the thermalization process in the self-interacting curvaton preheating scenario. We solve 
the evolution of the system with classical lattice simulations with a recently released symplectic 
PyCOOL program during the resonance and the early thermalization periods and compare the 
results to the inflaton preheating. After this we calculate the generated non-gaussianity with the 
AA'^ formalism and the separate universe approximation by running a large number of simulations 
with slightly different initial values. The results indicate a high level of non-gaussianity. We also use 
this paper to showcase the various post-processing functions included with the PyCOOL program 
that is available from https : //github . com/ j tksai/PyCDQL_ 
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I. INTRODUCTION 

The curvaton mechanism [IHT] is a much studied al- 
ternative to the standard inflationary paradigm for the 
origin of the observed primordial perturbations. The cur- 
vaton field is assumed to be light and subdominant dur- 
ing the inflation process and its contribution to the en- 
ergy density is significant only moments before its decay. 
This allows the inflation potential to have more natural 
properties [5j compared to the single field scenario while 
still leading to adiabatic perturbations consistent with 
the current observational data O [10] . 

Reheating of the universe is an important part of the 
early universe cosmology (for a review cf. |11|). In the 
curvaton scenario it is most often assumed that the cur- 
vaton field decays perturbatively into lighter degrees of 
freedom once the Hubble parameter is of the order of 
curvaton decay width F and thermalizes with the radia- 
tion that originates from the inflaton. It is however also 
possible that the universe reheated through a rapid and 
rather violent preheating process. This parametric res- 
onance was studied in ref. |12| in the curvaton scenario 
and the main conclusion was that in general it is quite 
similar to the preheating of the inflaton field. In ref. [13] 
it was further found that the curvaton resonance can lead 
to very high levels non-gaussianity. 

The curvaton potential in most of these studies is as- 
sumed to be of a quadratic type. As was noted in refs. 
[HI [T5] any deviations from this shape can lead to sig- 
nificant differences in the end results, especially in the 
level of generated non-gaussianity. Whereas these stud- 
ies were limited to the perturbative regime in the present 
paper we expand this analysis to the non-linear preheat- 
ing process. We limit the potential function of the curva- 
ton to the typical quadratic type with additional quartic 
self-interactions. We also assume that the curvaton field 
does not couple to other scalar fields in contrast to refs. 

[nma. 

We will study this self-interacting curvaton scenario 
with classical fields and lattice simulations from two dif- 



ferent perspectives. We will first concentrate on the ther- 
malization of the curvaton field during the resonance 
process. We will compare the results to the preheat- 
ing of inflaton that has been studied thoroughly in [151 - 
[5nj [2^^E5] with analytical and numerical methods. Af- 
ter this we will concentrate on the calculation of gener- 
ated non-gaussianity with the AA^ formalism [38]. This 
mainly numerical study will be done with the recently 
published symplectic PyCOOL program [37] (available 
from https://github.com/jtksai/PyCOOL). We also 
use this paper to showcase the numerous post-processing 
functions included with the program. 

This paper is organized as follows. In section II we 
present the curvaton model and the equations of motion. 
In section III we present the thermalization and non- 
gaussianity calculations and results. We conclude with a 
discussion in section IV. 



II. CURVATON MODEL 

We model the curvaton field with a simple polynomial 
potential function with quartic self-interactions 



K = ImW 



\Ka^ 



(1) 



where a is the curvaton field and A^. is the coupling con- 
stant of the curvaton self-interactions. Following [T3] we 
will set the initial energy density of the homogeneous ra- 
diation component equal to the potential energy of the 
inflaton 



v^ = \x^<l>' 



(2) 



where the coupling constant A^ is a free parameter and 
we set ~ rnpi, mpi being the reduced Planck mass. The 
curvaton fleld is effectively massless during inflation and 
hence it is require that 



3x^(7^ < h: 



(3) 
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where iJ* is the value of the Hubble parameter during 
inflation. 



After the inflation ends the curvaton field stays almost 
constant until it starts to oscillate around its minimum 
when the Hubble parameter has decreased close to the 
value of the effective mass of curvaton. In the usual per- 
turbative analysis the field would then start to decay into 
lighter particles once the Hubble parameter is roughly 
equal to the decay width of the curvaton. In this paper 
we are however more interested in the non-perturbative 
analysis meaning that the interaction terms in the poten- 
tial function (II]) now lead to the production of curvaton 
particles [3D] . The curvaton field is assumed to decay per- 
turbatively only long after the resonance period is over. 

The closely related reheating process of a self- 
interacting inflaton field has been studied previously in 
refs. [m [m HH 111130] of which the last two use a simi- 
lar interaction picture to this study. We will assume that 
the quartic term dominates the curvaton potential and 
hence initially we set a > m^j \f\^. In the opposite case 
the reheating process does not happen and the curvaton 
field does not thermalize. 

The creation of particles during this preheating has 
been studied extensively in [5D] in the case of massless 
inflation and we will cite the most relevant results here. 
The mode equation of the curvaton particles with wave 
number k can be written in terms of a more general Lame 
equation 



r 



cn2(,^,i=)jafc = 0, (4) 



which is valid for V\,, 



int — ^cr^X^ type interaction terms 
where x is another scalar field. This equation however re- 
duces to the mode equation of the quartic self-interaction 

2 

when |- = 3 [30 . We have here also defined a = acr, 
used prime to denote time derivative with respect to the 
scaled conformal time which is defined in terms of the 
physical time dt as df] = a^^\/\^aodt and cn(7y, -^) is 

the Jacobi cosine function. We have also used a rescaled 
wave number hP = k"^ /{X^gq), where the rescaled cur- 
vaton amplitude (Tq is measured at the end of inflation. 
The values of k^ and |- that will lead to production 
of particles can be read from the corresponding stabil- 
ity/instability chart that can be found for example in 
[3D]. 

We will now approximate the mode equation of the 
massive self-interacting curvaton particles with equation 

(4) with |- = 3 and we will also neglect the mass term 
wnich we assume to be small compared to the interaction 
term at least during the early part of the evolution. It is 
now easy to see from the stability/instability chart that 
the curvaton particles are produced at a band close to a 
rescaled momentum value of k^ ~ 1.6 which in terms of 
the comoving momentum reads 



l.eAcrCTo- 



(5) 



This resonant phase of particle production is followed 
by [19] a period of rescattering of the coherent curva- 
ton mode (fc — 0) and the created particles leading to a 
formation of multiple peaks in the spectrum of the field 
close to the harmonic frequencies of kp. After this the 
system enters a regime of turbulent dynamics [11 which 
is followed by a long period during which the field reaches 
the thermal state. 



A. Equations of motion 

We will solve the evolution of the system with a sym- 
plectic algorithm that is by design meant to conserve 
the energy of the system. Instead of solving the Euler- 
Lagrange equations of motion we will instead use the 
Hamiltonian equations that are split into explicitly in- 
tegrable pieces. Note that prime in the following equa- 
tions means derivative with respect to the conformal time 
dri = a~^dt. 

Starting from the Einstein-Hilbert action and after 
some simple Legendre transformations the Hamiltonian 
function of the system can be derived. Since we will solve 
the equations of motion numerically in a periodic comov- 
ing lattice the system needs to be discretized. We will 
use a second order accurate and fourth order isotropic 
stencils for the Laplacian operators derived in [28]. The 
discretized Hamiltonian function in conformal time then 
reads [27] 



n 
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2a^dx^ 

+ V{(l)i^S,--,4'N,x) , 



(6) 



where Pa is the canonical momentum of the scale factor 



a, Vl 



equals the size of the cubic lattice, dx the 



This is the only momentum band and the other suitable 
momentum values correspond to single points [3D]. 



spacing of the lattice mpi is the reduced Planck mass, 
TTi^s is the canonical momentum of field (t)i^s at position 
X — (xi, X2, Xa) in the lattice and D[(/)i^s]{x) is the Lapla- 
cian of field (pi^s at position x. Note also that the sum- 
mation is carried over all of the fields and all positions 
in the lattice. We have also incorporated homogeneous 
radiation p-y q a-^id non-relativistic matter Pm,o compo- 
nents into the system. It can be easily seen [57] that 
right hand side of the Hamiltonian ^ corresponds to 
the first Friedmann equation and is therefore conserved 
by the symplectic integrator. 

The Hamiltonian equations related to this Hamilto- 
nian now read for the scale parameter and its canonical 
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Similarly the equations of motion of scalar field i at grid 
point zread 



'^i.z- 



on 



a2 



= a 






t dV 



(8) 



which follow from equation ([6]) by differentiating under 
the summation sign and by summing over the coeffi- 
cients Cd(a) of the discretized Laplacian. When integrat- 
ing these equations we will first split them into explicitly 
integrable pieces and then use a suitable symplectic in- 
tegrator. 



III. NUMERICAL RESULTS 
A. Initial values 



these two cases. We have therefore used a comoving edge 
length 5/(12to) meaning that the comoving momenta are 
in the range 15.1m < k < 380 ttt- in the non-gaussianity 
results. 

The initial values for the curvaton field were chosen 
based on two criteria. In order for the quartic term to 
dominate in the potential function we simply set 



(To > 






(9) 



We also want the momentum band where the particle 
creation happens, i.e. Eq. ([s]), to be inside the lattice 
meaning that the parameters should be chosen such that 



n/tga; 



CTo 



(10) 



is neither too large nor too small. 

With these criteria in mind we used the following val- 
ues for the parameters: the mass of the curvaton is set 
to 1 X 10~^, initial curvaton field value ctq = 2 x 10"'*, 
curvaton self-interaction strength Acr = 1 x 10~^, initial 
radiation density A^ = 1 x 10^"'^^. The initial fractional 
energy density of the curvaton, fla,Q, corresponding to 
these values is of order ~ 10~^ (see Figure^. The mo- 
mentum band where the particle creation happens is ap- 
proximately at kp/m ^ 80. 



B. Output variables 



We use units where the reduced Planck mass ttipl is 
set to one. We will also use a general mass m — lO^^nipi 
to define the lattice, the initial radiation energy density 
and the time step dr]. The physical time is measured in 
units of m~^. We will use a conformal time step dr] = 
0.001/m in the simulations and solve the evolution until 
^phys m ~ 5000. 

The size of the lattice is limited by requirement that 
L < \/[aH) i.e. the comoving horizon is larger than 
the comoving lattice at all times. Otherwise the assump- 
tion that has been used when deriving equation (|6| that 
the metric is of the Friedmann-Robertson- Walker form 
—ds^ = a{'qY{—dif + dx^) would have to be adjusted to 
include also metric perturbations. 

We used two different lattice sizes to run the simu- 
lations: the thermalization study was done with 256'^ 
points whereas the non-gaussianity simulations were run 
on smaller 64"^ lattices that are roughly 46 times faster 
to solve. We set the comoving edge of the lattice to be 
5/(3r7T,) in the thermalization simulations meaning that 
the comoving momenta are in the range 3.8 rn < k < 
380 m which we calculate with the effective wave number 
fccff instead of the magnitude of the wave vector. In the 
non-gaussianity calculations with a smaller lattice size we 
are compelled to reduce either the infrared or the ultravi- 
olet resolution of the simulation. Simple numerical test 
runs show that the ultraviolet modes are more impor- 
tant for the evolution of the system to be consistent in 



Previous studies of the thermalization process after 
preheating have used a number of different variables to 
study and to illustrate the different phases of this process. 
The comoving number density and the related number 
density spectra are certainly some of the most interest- 
ing ones to use. There have been however a number of 
different definitions and ways to calculate these variables 
leading to slightly different results while the overall pic- 
ture of the thermalization process stays the same. In 
this study we use a definition for the number density Uk 
that was previously used in LATTICEEASY [55]. This 
is done by using conformal field values F^^c = o,Fk and 
conformal time to write the equations of motion of the 
Fourier modes of the fields in the form of a simple har- 
monic oscillator 



Pk,c + ^kPk,c 



0, 



(11) 



where -Ffc = L~^^'^Ff^ is the scaled Fourier mode of con- 
formal field af, L is the comoving length of the lattice 
and 



= fccff + a^™cff = fccff 






(12) 



is the comoving dispersion relation. Note that we have 
used the effective wave number k^g which is calculated 
from the discrete Fourier transform of the discretized 



Laplacian operator. The wave number is often how- 
ever calculated with the magnitude of the wave vector 
k'^ = fc^ + fcy + k^ as is done for example in LAT- 
TICEEASY. This method might however lead to inac- 
curate number density results [29, whereas the effective 
wave number takes properly the used discretization into 
account. We have also defined the effective mass rricff m 



equation (12 1 where the brackets denote an average over 



the lattice. The number density of the scalar particles 
can be then written in terms of the scaled modes F^ as 






(13) 



which is calculated by binning the data into spherical 
shells in the momentum space that are then averaged. 
We will also study the time evolution of the number of 
particles in the comoving lattice 



N{t) 



1 
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Ukd k, 



(14) 



which is calculated by summing over the non-averaged 
momentum bins. 

We are also interested in various energy density re- 
lated variables. We first define the energy density spectra 
based on the number density equation ([l3|) as 



Pk = ^knk 



(15) 



where now uJk = iok/o- is the physical dispersion relation. 
The energy density of a quanta at momentum k then 
simply reads e^ — k'^pk [23] ■ The energy and the pressure 
density of a scalar field in position space are defined as 
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(16) 
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We calculate the fractional energy densities from these 
expressions with 



a: = 



Pt 

Ptot 



(17) 



where ptot now includes all of the scalar fields and the ho- 
mogeneous radiation component. The equation of state 
is derived from ( 16 ) 



{P^)' 



(18) 



where the brackets denote averaging over the lattice. 

We are also interested in the statistical properties of 
the fields during the resonance process. In this study we 
will use the excess kurtosis which is defined as 



Pi „ 

72 = ^-3, 



(19) 



where p^ is the fourth moment about the mean and a is 
the standard deviation (not to be confused with the cur- 
vaton field). This quantity is mainly used to gauge how 



much the distribution of the curvaton field deviates from 
a gaussian one for which it is identically zero. A large 
value of kurtosis generally indicates that the distribution 
has more mass in the tails. 



C. Thermalization results 

We use a fourth order symplectic integrator to solve the 
evolution of the system in conformal time. The output is 
calculated after a constant number of integration steps. 
The moving averages presented in the figures are calcu- 
lated over these points meaning that when presented in 
physical time the length of the averaged period increases 
with time. We therefore use the term conformal moving 
average in the figures. 

The numerical accuracy during the simulation is shown 
in Figure [l] where we plot the absolute value of the resid- 
ual curvature 



K 



8TrG{p) 



3i/2 



-1 



(20) 



which we use to measure the conservation of Hamilto- 
nian([6]). As can be seen from the figure the algorithm is 
accurate to 10~^° level during the preheating phase. The 
error does increase with time but not substantially. 

The progress of the thermalization process is presented 
in Figure [2] where we plot the comoving number density 
as a function of time. As can be seen from the figure the 
number density initially stays close to a constant but as 
the resonance process starts the number density begins 
to increase exponentially. At iphys m ^ 40 the resonance 
ends and the system then enters the rescattering period. 
During this the number density reaches a short plateau 
phase after which it starts to gradually decrease mainly 
due to a lack of infrared resolution of the lattice. Overall 
the evolution of the number density is quite similar to 
the one witnessed in the chaotic inflation case jM]. 

Close inspection of the evolution of the number and en- 
ergy density spectra however tells a very different story 
when compared to the chaotic inflation. In the broad 
parametric resonance of the chaotic inflation the preheat- 
ing process is most efficient at creating particles with 
momentum values below a threshold value fc* |22j . In 
terms of the energy density of the quanta at momen- 
tum k the chaotic inflation potential usually leads to a 
spectrum with one peak at the infiaton particle energy 
spectrum that broadens with time and shifts to higher 
comoving momentum values with time [53j. In the self- 
interacting curvaton case the particle creation happens 
initially at the resonance band calculated in eq. ([5]) which 
can be seen in Figure |3] as a formation of a clear peak at 
k/m ~ 80. This phase is however followed shortly by ex- 
citation of curvaton particles at a series of different bands 
indicating that the system has entered the rescattering 
period [19]. Note that this part of the process is quite 
sensitive to the initial values: at larger initial radiation 
densities or smaller curvaton self-interaction values it is 
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FIG. 1. The evolution of the numerical error during the sim- 
ulation (thin line) and its conformal moving average (thick 
line). Notice that the averages are calculated over the values 
at different output points which are written after a constant 
number of integration steps in conformal time. Therefore the 
physical time over which the moving average is calculated 
varies. Notice also that the large gaps at fphys/"i ~ 1400 and 
tphys/m ~ 3800 in the graph are an artifact of this used sam- 
pling. Please see the online version of this article for color 
figures. 
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FIG. 2. The evolution of the comoving number density of 
the curvaton particles. Different phases of the process are 
clearly distinguishable: exponential resonance period at 1 < 
tphys/m < 40, the rescattering phase 40 < tphys/m < 200 and 
the final turbulence period. 
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FIG. 3. The evolution of the number density spectrum. 
Note that in the figure the color evolves with time and the 
red curves are calculated close to the end of the simula- 
tion whereas the blue ones (at the bottom) are evaluated 
at tphys m = 0. The time difference between the spectra is 
roughly tphys m, ~ 20. We have also included a power law 
fit Uk ~ k~^ of the final spectrum with p ~ 3/2 as a black 
dashed curve in the figure. Please see the online version of 
this article for color figures. 




FIG. 4. The evolution of the energy density spectrum of the 
quanta at momentum k. Note that in the figure the color 
evolves with time and the red curves are calculated close to the 
end of the simulation whereas the blue ones (at the bottom) 
are evaluated at iphys m — 0. The time difference between the 
spectra is roughly iphys m ~ 20. Please see the online version 
of this article for color figures. 
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possible to stop this process before the other peaks start 
to form. 

The shape of the number density spectrum at the end 
of the simulation is visible in FigurelSlas a red curve. The 
observed peaks have leveled out except for small residual 
hills. Other notable feature is that the spectrum is ele- 
vated at smaller momentum values. This final shape also 
appears to be quite stable in the sense that it does not 
change considerably during the last stages of the simula- 
tion. To compare this to a thermal boson spectrum we 
have fitted the data to the usual Rayleigh- Jeans approx- 
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(21) 



where T is the comoving temperature of the boson field in 
thermal equilibrium and /i is the corresponding chemical 
potential. The best least squares fit (not shown in the fig- 
ure) results in /i ~ a TTioff and T ^ a mpi which strongly 
indicates that the system is non-thermal. A power law 
function n^ ex fc~^ with p ~ 3/2 seems to follow the shape 
of the spectrum more closely until an exponential cut-off 
at high momentum values. Similar result was previously 



presented in the case of self-interacting massless inflaton 
field in ref. [21] where the evolution of the spectra during 
the turbulence period was in addition found to be self- 
similar. Although we were unable to verify this with the 
curvaton model the results indicate that the curvaton is 
not at thermal equilibrium at the end of the simulation. 
Assuming that the eventual thermalization of the curva- 
ton happens through the quartic interactions and that 
it is not coupled to other fields the corresponding decay 
rate reads 
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(22) 



which leads to a rather low reheating temperature of a 
few MeV. 

Another perspective to the resonance process can be 
seen in Figure |4] where we plot the evolution of the energy 
density of the quanta at momentum k with the quantity 
k^Pk- As is evident from the graph most of the curvaton 
particles are created at five different harmonic momen- 
tum bands. As time evolves the series of peaks smoothen 
as the thermalization process progresses and the energy 
density of the particles propagates toward higher momen- 
tum values. The final state in this case is very different 
from the one seen in the chaotic inflation 1231. 
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FIG. 5. The evolution of the equation of the state of the 
curvaton Ua i.e. Eq. (181 during the simulation. Note that 



the thin line is the variable and the thick line is the conformal 
moving average. Notice also that the large gaps at tphys/m ~ 
1400 and tpi^ys/m ~ 3800 in the graph are an artifact of the 
output sampling. 

We are also interested in the evolution of the equation 
of the state of the curvaton during the thermalization 
process. As can be seen from Figure [5] the system is 
initially highly relativistic and oscillates rapidly. This 
oscillatory phase corresponds to the exponential increase 
in the comoving particle number density seen in Figure 
[2j As the system evolves the equation of state starts to 
decrease but at the end of simulation its average is still 
close to a value of 0.1 indicating that the system is not 
yet non-relativistic. 

In Figure [6] we plot the comoving effective mass a rrics 
in units of m which we calculate as an average over the 
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FIG. 6. The evolution of the comoving effective mass a rrics 
during the simulation. Note that the thin line is the variable 
and the thick line is the conformal moving average. 
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FIG. 7. The evolution of the fraction of relativistic curvaton 
particles, i.e. for which k > a rricft, during the simulation. 
Note that the thin line is the variable and the thick line is the 
conformal moving average. 



lattice. The early stages are in this case also highly oscil- 
latory which is followed by a period of gradual increase 
due to the expansion of the universe. During the reso- 
nance and the rescattering periods the comoving effective 
mass stays almost constant and it starts to grow only af- 
ter the mass term starts to dominate at iphys 'tt- ^ 200. 

In Figure [7| we show the fraction of curvaton particles 
that are relativistic i.e. for which k > a nies and the 
homogeneous mode is not included in the calculations. 
The figure shows that the created curvaton particles are 
highly relativistic during the simulation with a final value 
close 65 percent. Notice that the discrepancy between 
Figures [5] and [T] is caused by the coherent curvaton field 
that still gives a significant contribution to the energy 
and pressure densities of the curvaton component at the 
end of the simulation. 

Yet another aspect of the evolution of the curvaton is 
seen in Figure [8] where we plot the fractional energy den- 
sity of the curvaton during the simulation. Initially it 
evolves in tandem with the homogeneous radiation com- 
ponent up to time tm ~ 10 after which its fraction of 
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FIG. 8. The evolution of the fractional energy density of the 
curvaton fla as a function of time. 



energy density starts to grow steadily as its equation of 
state starts to approach that of matter. 
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FIG. 9. The evolution of excess kurtosis i.e. Eq. (19 1 during 



the simulation. Note that time is given in logarithmic units. 

We have also included a plot of the excess kurtosis 72 
during the thermalization process in Figure [91 In the 
early highly oscillatory preheating phase the system is 
also highly non-gaussian. However as the simulation pro- 
gresses the curvaton field starts to return to gaussian. 
This behavior is very similar to the one observed in the 
parametric resonance of the chaotic inflaton field [24] . 
The skewness of the curvaton field shows a very simi- 
lar trend and we have therefore omitted the graph of its 
evolution. 



D. Non-gaussianity calculations 

The possible generation of non-gaussianity during the 
curvaton thermalization process is an interesting and a 
timely question in cosmology [lUl I5TH55] . To calculate 
this we will use the A7V formalism based on the separate 
universe approach [38J that has been previously applied 
successfully to different parametric resonance scenarios 
[T51 [531 - I5S] . In the separate universe approach different 



patches of the universe that are separated by more than 
a Hubble distance are presumed to evolve independently 
of each other. Assuming also that each Hubble volume 
is isotropic and homogeneous they can be approximated 
to be separate Friedmann- Robertson- Walker 'universes'. 
The evolution of these patches is solved with the lattice 
simulation method as in the previous section. 

The curvature perturbation on scales larger than the 
Hubble horizon is defined as 



C = (51na| 



H, 



(23) 



where the difference in the scale factor is calculated at 
a hypersurface of constant Hubble parameter H. The 
scale factor is normalized to be one at the start of the 
curvaton thermalization process. We will vary the ho- 
mogeneous value of the curvaton field with superhorizon 
fiuctuations from one patch to another. This will cause 
slight variations in the value of the curvature perturba- 
tion C. For small perturbations Sa equation (23) is often 
expanded as 
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where the primes are derivatives calculated with respect 
to the curvaton value at the end of inflation on hypersur- 
faces of constant Hubble parameter H. The spectrum of 
the curvature perturbation can be written with this as 
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where Po- is the spectrum of the curvaton field. Following 
[131 133 we will use 
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which is valid for massless fields during infiation. Nk{~ 
60) here measures how many number of e-foldings before 
the end of inflation mode k = aHk left the Hubble hori- 
zon. The local non-gaussinity parameter can be defined 
also in terms of the coefficients of equation ( 24 ) [37] as 
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To calculate the non-gaussianity in the curvaton sce- 
nario we will apply the method presented in (13) with 
minor modifications. We will write the energy density as 
a combination of the relativistic radiation and the curva- 
ton component which we assume to behave like matter 



Prof 



^..ef( — )' + (l-r,ef)(— )1, (28) 



where the fractional energy density of curvaton r^ef = 
^cr.rof, scale factor aref and energy density pref are cal- 
culated at a reference point defined after the resonance 
period of the curvaton. 



We will assume that the curvaton stays subdominant 
during its evolution and decays perturbatively when the 
Hubble parameter H is of the order of the decay width F. 
We will use the sudden decay approximation by assuming 
that this decay is instantaneous. The value of the decay 
width is unknown meaning that the energy density pdocay 
and the fractional energy density rjecay at the moment of 
decay are free parameters limited by observational data, 
namely the amplitude of the curvature perturbations. By 



now taking logarithms on both side of equation ( 28 ) , ex- 
panding the right side in series with respect to r,.cf and 
rearranging the terms the logarithm of the scale factor 
reads 



(29) 
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where 
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(30) 



and C = r/r-rci— 1. The curvature perturbation can now 
be written as 



C((To) = Ina(ao) — lna((To) = 5hia^ci + CSr^:, 



(31) 



where we have written explicitly the dependence on the 
curvaton value at the end of inflation. We have also ne- 



glected the energy density terms from equation ( 29 1 since 



the calculations are done on a constant H hypcrsurface 
on which also the energy density is constant by the Fried- 
mann equations. 

We will now assume that the logarithm of the scale 
factor and the fractional energy density can be expanded 
in terms of the superhorizon fluctuations of the homoge- 
neous curvaton values similarly to equation (24): 



(32) 



Inaref(o-o) =lnarof(a-o) + lna;.j,f((To - ctq) 

»^rof(o-o) =?'rot(CTo) + ?'rcf('^0 " C^o) 

+ 2^"of('^0 -C^o)^ 

where (Tq = cro+<^fo a-nd Suq is a superhorizon fluctuation 
of the initial curvaton value. Equations (32 1 are fitted 



to the simulation data to get numerical values for the 



polynomial coefficients InaJ, 



rcf' 



ln<of> <of 



and 



rcf' 



For 



the amplitude of the curvature perturbation spectrum 
( [25| ) to be consistent with the WMAP observations [9], 
Pq ~ 2.4 X 10~^, the unknown fractional energy density 
of the the curvaton at the moment of decay can be solved 
|13j in terms of the power spectrum amplitudes and the 
polynomial coefficients 



^dccay ^rcf ^ ^ 




ln( 
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(33) 



The non-gaussianity parameter (27) can be written sim- 
ilarly [13] as 







ln<cf ■ (34) 
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FIG. 10. The fractional energy density r^^t calculated at the 
reference value of the Hubble parameter Href as a function of 
the initial homogeneous value of the curvaton. The contin- 
uous line in the graph represents a second order polynomial 
least squares fit to the data. We have also included standard 
error limits of the simulation data. 
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FIG. 11. The difference of the logarithm of the scale factor 
a calculated at il/rof as a function of the initial homogeneous 
value of the curvaton. The continuous hue in the graph repre- 
sents a second order polynomial least squares fit to the data. 
We have also included standard error limits of the simulation 
data. 

The Monte Carlo simulations were run with the initial 
values that were used in the thermalization analysis. As 
mentioned previously we used a smaller lattice size of 64^ 
points in order to shorten the overall simulation runtime 
drastically (roughly 46 times faster). For the reference 
point where the different quantities are calculated we use 
i^rcf = 4 X lO^^JTi which in terms of physical time cor- 
responds to tm ~ 1260. The actual value is determined 
by interpolating around TJrcf • The range of homogeneous 
curvaton values over which the simulations need to be 
run is determined by the variance of the curvaton values 



at the end of inflation. For inflation potential (pi) and 
curvaton spectrum ( 26 ) this reads [13] 



{Sa' 



9tt 



2^™PL^^0 



(35) 



where iVo ~ 60 is the number of e-foldings after the 
largest currently observable scales left the horizon. The 
range of curvaton initial values then reads 



O'O - o'^CTo < o-Q < ctq + i:,Scro, 



(36) 



where Sao = \/{S^ ~ 9.9 x 10"^ and cto = 0.0002. We 
take 41 equidistant points from this range and use as the 
initial homogeneous curvaton values. At each point the 
simulations are solved with different random field per- 
turbations 35 times to get the necessary statistics. Note 
that these subhorizon perturbations are generated with 
a convolution based algorithm presented in [5S] . The to- 
tal simulation runtime with these selections is roughly 25 
hours when using a Nvidia Tcsla C2050 computing card. 




-0.00005 



that fla ~ Pa/ P-y ~ a during this period the value of 
the Hubble parameter at decay can be calculated to be 
roughly 0.1 eV which translates to a reheating temper- 
ature Tich ^ 1 GeV which is considerably higher than 
the result of the previous section. The non-gaussianity 
variable can be calculated from equation (34) or by fit- 
ting C directly with equation (24 1. The results are 
/nl = 2980 ± 644 and /nl = 2976 ± 1481 respectively at 
95% confidence level. When compared to the results of 
a two field curvaton resonance model |13j the quadratic 
polynomials follow more closely the general trend of the 
data. Despite this the calculated level of non-gaussianity 
is still very high and the current observational limit ^ID] 
—9 < /nl = 111 at 95% confidence level rules out the 
model with the current parameter values. 

This large level of non-gaussianity is mainly caused by 
the magnitude of the second order coefficient r'.'^j and 
the smallness of the first order coefficient r'^^^ in equa- 
tion (34 1. An easy remedy to this would be to use a 



smaller curvaton self-interaction strength which would 
lead to a more linear evolution of the fractional energy 



density of the curvaton in Figure 10 This might however 
cause some thermalization related problems mentioned 
briefly in the previous section: for smaller values of self- 
interaction coupling strength the rescattering phase after 
the resonance period was found to be very weak and lim- 
ited and the final shape of the number spectrum exhibit 
a clear peak at fc ^ fcp. The created comoving number 
density of the particles in this case would be also orders 
of magnitude smaller than with the current values. 



IV. DISCUSSION AND CONCLUSIONS 



FIG. 12. The curvature perturbation (" as a function of the 
initial homogeneous value of the curvaton. The continuous 
line in the graph represents a second order polynomial least 
squares fit to the data. We have also included standard error 
limits of the simulation data. 



The main results of the simulations are presented in 
Figures [T0p2| and in Table [H In Figures [lU] and [Tl] we 
have the fractional energy density r^cf and the difference 
of the logarithm of the scale factor a calculated at the ref- 
erence value of the Hubble parameter TJref as a function 
of the initial homogeneous value of the curvaton. The 



curvature perturbation C, calculated with formula (31) is 
given in Figure [T2| We have also included least square 



fits of the equations (32) in the graphs with the corre- 
sponding polynomial coefficients given in Table |lj Note 
that we have also included the confidence intervals of the 
parameters at 95 % level which were derived from the 
fitting results given by Mathematica. 

With these results the curvaton fraction at decay reads 
''decay = 0.037 ± 0.0024. Assuming that the radia- 
tion stays dominant after the end of the simulation and 



We have studied the self-interacting curvaton scenario 
with classical fields and lattice simulations from two dif- 
ferent perspectives. First we concentrated on the ther- 
malization process during and after the preheating phase. 
The results indicate that in the current curvaton scenario 
the overall evolution of system follows closely the previ- 
ously studied self-interacting inflaton model. We found 
that during the resonance period curvaton particles were 
created at a predicted resonance band and in the ensu- 
ing rescattering phase the spectrum developed peaks at 
harmonic frequencies related to the momentum values of 
the resonance band. The final state of the curvaton field 
could be characterized as a pre-thermalized one. 

After this we concentrated on the calculation of the 
generated non-gaussianity during the resonance. We em- 
ployed and adapted a previously presented method T3J to 
the self-interacting curvaton scenario. When compared 
to the broad resonance of curvaton jT3] the simulation 
data was found to be a better fit to the used quadratic 
approximation of the curvature perturbation. The used 
parameter values were however rule out by the current 
observational limits and were found to be unphysical. 
There might however be regions in the parameter space 
that could lead to non-gaussianities consistent with the 





1 


(To — (To 


((To - aof 


lnaref(o"o) — Inarcf('To) 


(6.465 ±5.058) x 10"'° 


(1.793 ±0.1155) X lQ-2 


7885 ± 4430 


rref((To) 


(2.693 ± 0.002022) x lO"*^ 


(7.170 ±0.4620) X 10"^ 


31530 ± 17720 


c 


(8.868 ± 6.947) x 10-« 


246.3 ±15.87 


(1.083 ±0.6086) X 10* 



TABLE I. List of least squares fit results for second order polynomials lnarcf((To), ''rof(iTo) and ^ i.e. equations (32 I and (241 
respectively. In the columns we have given the coefficients of different powers of (To — (To. We have also given the confidence 
intervals for the parameters at 95 % level. 



observations. This would however take more comput- 
ing resources that were available while doing this paper. 
One option would be to make a distributed version of 
PyCOOL that would systematically scan the parameter 
space for suitable initial values. Another interesting pos- 
sibility would be to study the generation of gravitational 
waves during the curvaton resonance which would likely 
give additional limitations on the curvaton model. This 
could be done easily with a recently updated version of 
PyCOOL and we leave it for future work. 
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